home *** CD-ROM | disk | FTP | other *** search
/ Travel to Space / Travel to Space.iso / dos_prog / astronmy / geosat / geosat.pas < prev   
Pascal/Delphi Source File  |  1990-10-16  |  14KB  |  461 lines
  1. {$A+,B-,D+,E+,F-,I-,L+,N+,O-,R-,S-,V-}
  2. {$M 16384,0,655360}
  3. (*
  4.                              Program 'GEOSAT'
  5.                Adapted to Turbo Pascal on October 13, 1990
  6.                              by Philip Miller
  7.                                    from
  8.                          the Quick Basic Program
  9.                     Copyright (C) 1982 by David Eagle
  10.              Public domain for the IBM-PC on October 8, 1986
  11.  
  12.          Computes direction of geosynchronous satellites
  13.          Azimuth angle   : Positive clockwise from north ( degrees )
  14.          Elevation angle : Positive above the horizon ( degrees )
  15.  
  16. *)
  17.    { ********************************************************** }
  18.  
  19. Program GeoSat;
  20.  
  21. Uses OpCrt,
  22.      OpString,
  23.      OpDate;
  24.  
  25. {$I FLOAT.INC}        (* determines float type real or double *)
  26.  
  27. Type
  28.   SatRecType = Record
  29.     Bird    : String[12];
  30.     Code    : String[4];
  31.     Geolong : String[5];
  32.     Az      : String[13];
  33.     Alt     : String[13];
  34.   end;
  35.  
  36. VAR
  37.  A,B,C,D,R : FLOAT;
  38.  I, J, K   : Integer;
  39.  
  40.  HalfPi    : Float;
  41.  DTR       : Float;
  42.  RTD       : Float;
  43.  
  44.  LV        : ARRAY[1..3] OF FLOAT;
  45.  UP        : ARRAY[1..3] OF FLOAT;
  46.  
  47.  OBSLAT    : FLOAT;
  48.  OBSLONG   : FLOAT;
  49.  OBSALT    : FLOAT;
  50.  XOBS      : FLOAT;
  51.  YOBS      : FLOAT;
  52.  ZOBS      : FLOAT;
  53.  AZIMUTH   : FLOAT;
  54.  ELEVATION : FLOAT;
  55.  SATLONG   : Float;
  56.  
  57.  OBSLATstr    : String;
  58.  OBSLONGstr   : String;
  59.  OBSALTstr    : String;
  60.  AZIMUTHstr   : String;
  61.  ELEVATIONstr : String;
  62.  SATLONGstr   : String;
  63.  
  64.  SatRec       : Array[ 1..30 ] of SatRecType;
  65.  MaxRecs      : Byte;
  66.  
  67.  HelpAttr : Word;
  68.  HeadAttr : Word;
  69.  TitleAttr: Word;
  70.  DispAttr : Word;
  71.  LabelAttr: Word;
  72.  
  73. CONST
  74.  RREQ   = 1 / 637814;
  75.  FLAT   = 1 / 298.257;
  76.  HGEO   = 6.6107426;
  77.  SATFILE = 'GEOSAT.DAT';
  78.  
  79.    { *************************************************************** }
  80.  
  81.   Procedure SetColors;          { This a TURBO-POWER(r) library routine }
  82.     {-Choose attributes for display}
  83.     Begin
  84.       case CurrentMode of
  85.         3 :
  86.           Begin
  87.             HelpAttr := $1A;   { all regular screen colors }
  88.             HeadAttr := $1F;   { top header colors }
  89.             TitleAttr:= $1B;   { scroll bar colors }
  90.             DispAttr := $1E;   { category header colors }
  91.             LabelAttr:= $1C;
  92.             TextAttr := $1E;   { large frame color }
  93.           end;
  94.         ELSE
  95.           Begin
  96.             HelpAttr := $07;
  97.             HeadAttr := $0F;
  98.             TitleAttr:= $0F;
  99.             DispAttr := $07;
  100.             LabelAttr:= $0F;
  101.             TextAttr := $07;
  102.           end;
  103.       end;
  104.     end;
  105.  
  106.  
  107.  Procedure Init;
  108.    Begin
  109.      SetColors;
  110.      HalfPi := 0.5 * PI;
  111.      DTR    := PI / 180;
  112.      RTD    := 180 / PI;
  113.      FillChar( SatRec, Sizeof( SatRec ), #0 );
  114.      A := 0;
  115.      B := 0;
  116.      C := 0;
  117.      D := 0;
  118.      R := 0;
  119.    end;
  120.  
  121.  
  122.  Function Sgn( ww : float ) : float;
  123.    Begin
  124.      IF WW < 0 THEN
  125.        SGN := -1
  126.       ELSE
  127.        SGN := +1;
  128.    end;
  129.  
  130.  Procedure ATAN3( AA : float;
  131.                   BB : float;
  132.               var CC : float ); { 4-quadrant arc tangent subroutine }
  133.  
  134.    Begin
  135.      IF ABS(AA) < 1E-8 THEN
  136.        Begin
  137.          CC := (1-SGN(BB)) * HalfPi;
  138.          EXIT;
  139.        end
  140.      ELSE
  141.        CC := ( 2 - SGN(AA)) * HalfPi;
  142.  
  143.      IF ABS(BB) < 1E-8 THEN
  144.        EXIT
  145.      ELSE
  146.        CC := CC + SGN(AA) * SGN(BB) * ( ABS( ArcTan( AA/BB ) ) - HalfPi );
  147.    end;
  148.  
  149.  Function AzQuad16( vAZ : float ) : string;
  150.    Begin
  151.      IF (vAZ >= 0 ) and (vAZ < 11.25 )   OR
  152.         (vAZ >= 371.25) and (vAZ <= 360) THEN
  153.        AzQuad16 := ' (N)  '
  154.      ELSE IF ( vAZ >=  11.25) and (vAZ <  33.75) THEN
  155.        AzQuad16 := ' (NNE)'
  156.      ELSE IF ( vAZ >=  33.75) and (vAZ <  56.25) THEN
  157.        AzQuad16 := ' (NE) '
  158.      ELSE IF ( vAZ >=  56.25) and (vAZ <  78.75) THEN
  159.        AzQuad16 := ' (ENE)'
  160.      ELSE IF ( vAZ >=  78.75) and (vAZ < 101.25) THEN
  161.        AzQuad16 := ' (E)  '
  162.      ELSE IF ( vAZ >= 101.25) and (vAZ < 123.75) THEN
  163.        AzQuad16 := ' (ESE)'
  164.      ELSE IF ( vAZ >= 123.75) and (vAZ < 146.25) THEN
  165.        AzQuad16 := ' (SE) '
  166.      ELSE IF ( vAZ >= 146.25) and (vAZ < 168.75) THEN
  167.        AzQuad16 := ' (SSE)'
  168.      ELSE IF ( vAZ >= 168.75) and (vAZ < 191.25) THEN
  169.        AzQuad16 := ' (S)  '
  170.      ELSE IF ( vAZ >= 191.25) and (vAZ < 213.75) THEN
  171.        AzQuad16 := ' (SSW)'
  172.      ELSE IF ( vAZ >= 213.75) and (vAZ < 236.25) THEN
  173.        AzQuad16 := ' (SW) '
  174.      ELSE IF ( vAZ >= 236.25) and (vAZ < 258.75) THEN
  175.        AzQuad16 := ' (WSW)'
  176.      ELSE IF ( vAZ >= 258.75) and (vAZ < 281.25) THEN
  177.        AzQuad16 := ' (W)  '
  178.      ELSE IF ( vAZ >= 281.25) and (vAZ < 303.75) THEN
  179.        AzQuad16 := ' (WNW)'
  180.      ELSE IF ( vAZ >= 303.75) and (vAZ < 326.25) THEN
  181.        AzQuad16 := ' (NW) '
  182.      ELSE IF ( vAZ >= 326.25) and (vAZ < 348.75) THEN
  183.        AzQuad16 := ' (NNW)';
  184.    end;
  185.  
  186.  
  187.  
  188. Procedure OBSVECTOR;  { Observer position vector subroutine }
  189.  
  190.  VAR
  191.   E : Float;
  192.   F : Float;
  193.   G : Float;
  194.   H : Float;
  195.  
  196.  Begin
  197.   A := OBSLAT * DTR;
  198.   B := SIN(A);
  199.   C := COS(A);
  200.   D := 1 - FLAT;
  201.   E := SQRT( 1 - (2 * FLAT - SQR(FLAT) ) * SQR(B) );
  202.   F := 1 / E + OBSALT * RREQ;
  203.   G := SQR(D)/E + OBSALT * RREQ;
  204.   H := -OBSLONG*DTR;
  205.   XOBS := F * C * COS(H);
  206.   YOBS := F * C * SIN(H);
  207.   ZOBS := G*B;
  208.  end;
  209.  
  210.  
  211.  Procedure TOPOCENTRIC;         { Convert to topocentric coordinates }
  212.    VAR
  213.      MATRIX    : ARRAY[1..3,1..3] OF FLOAT;
  214.  
  215.    Begin
  216.      A := SIN( OBSLAT*DTR );
  217.      B := COS( OBSLAT*DTR );
  218.      C := SIN( -OBSLONG*DTR );
  219.      D := COS( -OBSLONG*DTR );
  220.      MATRIX[1,1] := A*D;
  221.      MATRIX[1,2] := A*C;
  222.      MATRIX[1,3] := -B;
  223.      MATRIX[2,1] := -C;
  224.      MATRIX[2,2] := D;
  225.      MATRIX[2,3] := 0;
  226.      MATRIX[3,1] := B*D;
  227.      MATRIX[3,2] := B*C;
  228.      MATRIX[3,3] := A;
  229.  
  230.      FOR I := 1 TO 3 Do
  231.        Begin
  232.          A := 0;
  233.          FOR J := 1 TO 3 Do
  234.            A := A + MATRIX[I,J] * UP[J];
  235.          LV[I] := A;
  236.        end;
  237.  
  238.      IF ABS( LV[3] ) > 1 THEN
  239.        ELEVATION := SGN( LV[3] ) * 90
  240.      ELSE
  241.        ELEVATION := ( ArcTan( LV[3] / SQRT( 1 - SQR(LV[3]) ) ) ) * RTD;
  242.  
  243.      A :=  LV[2];
  244.      B := -LV[1];
  245.      ATAN3(A,B,C);
  246.      AZIMUTH := RTD*C;
  247.  
  248.      ELEVATIONstr := Form( '###.##', ELEVATION );
  249.      AZIMUTHstr   := Form( '###.##', AZIMUTH );
  250.    end;
  251.  
  252.  Procedure PAUSE;           { Screen pause }
  253.   var
  254.    Dummy : char;
  255.  
  256.    Begin
  257.      gotoXY( 25, 24 );
  258.      Dummy := #0;
  259.      FastCenter( '< Press Any Key To Continue... >', 25, HeadAttr );
  260.      Dummy := ReadKey;
  261.      Repeat Until dummy <> '';
  262.    end;
  263.  
  264.  
  265.  Procedure INTRODUCTION;    { Program introduction subroutine }
  266.  
  267.    Begin
  268.      ClrScr;         { Uses some TURBO-POWER(r) ROUTINES }
  269.      FastWrite( 'Program "GEOSAT" is an interactive TURBO-PASCAL program which can      ',  4,5,HelpAttr );
  270.      FastWrite( 'be used to determine both the azimuth and elevation siting angles      ',  5,5,HelpAttr );
  271.      FastWrite( 'for a satellite dish located in North America pointing at the various  ',  6,5,HelpAttr );
  272.      FastWrite( 'geosynchronous satellites. The program takes into account the          ',  7,5,HelpAttr );
  273.      FastWrite( 'effect of the flattening or ''oblateness'' of the earth on the         ',  8,5,HelpAttr );
  274.      FastWrite( 'observer location on the earth. It also compensates for ''non-sea''    ',  9,5,HelpAttr );
  275.      FastWrite( 'level observation and dish sites. These features provide the user      ', 10,5,HelpAttr );
  276.      FastWrite( 'with very accurate pointing angles.                                    ', 11,5,HelpAttr );
  277.  
  278.      FastCenter( 'USER INPUTS AND SELECTIONS'                                            ,13,  HeadAttr );
  279.      FastWrite( 'The program GEOSAT will prompt the user for geographic coordinates to ' ,15,5,HelpAttr );
  280.      FastWrite( 'convert geosyncronous longitude to topocentric (earth) coordinates    ' ,16,5,HelpAttr );
  281.      FastWrite( 'The following describes the necessary input.                          ' ,17,5,HelpAttr );
  282.      PAUSE;
  283.  
  284.      ClrScr;
  285. {}   FastWrite( 'Observer Latitude (degrees):'                                         , 4,5,HeadAttr );
  286.      FastWrite( 'The  user  should respond with the geographic latitude in decimal  '  , 5,5,HelpAttr );
  287.      FastWrite( 'degrees taking note of the sign convention. For example, an        '  , 6,5,HelpAttr );
  288.      FastWrite( 'observer at 36 degrees, 30 minutes north latitude would input 36.5.'  , 7,5,HelpAttr );
  289.  
  290. {}   FastWrite( 'Observer West Longitude (degrees):'                                   , 9,5,HeadAttr );
  291.      FastWrite( 'Input the observer West longitude in decimal degrees. This will    '  ,10,5,HelpAttr );
  292.      FastWrite( 'always be a positive number.                                       '  ,11,5,HelpAttr );
  293.  
  294. {}   FastWrite( 'Observer altitude (meters):'                                          ,14,5,HeadAttr );
  295.      FastWrite( 'Input your altitude (height from sea level) at the above location  '  ,15,5,HelpAttr );
  296.      FastWrite( 'in meters. Sites above sea level are input as positive numbers and '  ,16,5,HelpAttr );
  297.      FastWrite( 'those below sea level are negative numbers.                        '  ,17,5,HelpAttr );
  298.      PAUSE;
  299.  
  300.      Clrscr;
  301. {}   FastCenter( 'PROGRAM OUTPUT'                                                    , 4,  HeadAttr );
  302.      FastWrite( '"GEOSAT" will output the observer-centered or topocentric azimuth  ',  6,5,HelpAttr );
  303.      FastWrite( 'and elevation angles at the Satellite Dish location for all the    ', 7,5,HelpAttr );
  304.      FastWrite( 'North American Satellites. Azimuth is measured positive clockwise  ', 8,5,HelpAttr );
  305.      FastWrite( 'from North and elevation is positive above the observer''s horizon.', 9,5,HelpAttr );
  306.      FastWrite( 'For example, an Azimuth reading of 90° is EAST and 180° is SOUTH.  ',10,5,HelpAttr );
  307.      FastWrite( 'Elevation angles which are negative indicate that the satellites   ',11,5,HelpAttr );
  308.      FastWrite( 'are below the horizon and therefore cannot be seen from an         ',12,5,HelpAttr );
  309.      FastWrite( 'observer''s location. Angles are printed in decimal degrees.       ',13,5,HelpAttr );
  310.      PAUSE;
  311.    end;
  312.  
  313.  
  314.  Procedure Calculate_Obs_Position;     { Calculate observer's position}
  315.   var
  316.    X3   : Float;
  317.    Y3   : Float;
  318.    Z3   : Float;
  319.    XSAT : Float;
  320.    YSAT : Float;
  321.  
  322.    Begin
  323.     OBSVECTOR;  { Calculate Satellite's Position }
  324.  
  325.     A    := -SATLONG*DTR;
  326.     XSAT := HGEO * COS(A);
  327.     YSAT := HGEO * SIN(A);
  328.  
  329.                 { Calculate Pointing Angles to Satellite }
  330.  
  331.     X3 := XSAT - XOBS;
  332.     Y3 := YSAT - YOBS;
  333.     Z3 := -ZOBS;
  334.  
  335.     R := SQRT( SQR(X3) + SQR(Y3) + SQR(Z3) );
  336.  
  337.     UP[1] := X3/R;
  338.     UP[2] := Y3/R;
  339.     UP[3] := Z3/R;
  340.  
  341.     TOPOCENTRIC;
  342.   end;
  343.  
  344.  Procedure Display_Results;
  345.    Begin                        { Print or display results}
  346.      ClrScr;
  347.      FastWrite( 'Observers Location:'         , 7, 61, TitleAttr );
  348.      FastWrite( 'Latitude:  ' + OBSLATstr     , 8, 61, HeadAttr );
  349.      FastWrite( 'West Long: ' + OBSLONGstr    , 9, 61, HeadAttr );
  350.      FastWrite( 'Altitude:  ' + OBSALTstr     ,10, 61, TitleAttr );
  351.      FastWrite( 'SAT NAME    CODE    GEOSYNCH    AZIMUTH         ELEVATION', 1,4, LabelAttr );
  352.  
  353.      FOR J := 1 TO MaxRecs DO        { was 20 - now determined by reading }
  354.        With SatRec[J] DO
  355.          FastWrite( Bird + '  ' + Code  + '  ' + Geolong + '°      '
  356.                          + Az   + '   '+ Alt  , J+1, 4, TextAttr );
  357.      PAUSE;
  358.    end; {Display_Results}
  359.  
  360.  
  361. Procedure MainProgram;
  362.   Label
  363.    Location,
  364.    Satellite,
  365.    Main;
  366.  
  367.   Var
  368.    INTRO : Char;
  369.    REPLY : Char;
  370.    FP    : Text;
  371.  
  372.   Begin
  373.     Init;     { Init variables, set colors }
  374.  
  375.     ClrScr;
  376.     FastCenter( ' GEO-SAT -- SYNCRONOUS SATELLITE POSTIONS ', 2,HeadAttr );
  377.     FastCenter( ' (C) Copyright 1990 by Philip Miller ', 3,TitleAttr );
  378.     FastWrite( 'Program introduction ( <Y>es, <N>o )', 5,5,TitleAttr );
  379.     gotoXY( 45, 5 );
  380.     INTRO := ReadKey;
  381.     IF INTRO in ['y','Y'] THEN
  382.      Begin
  383.        HIDDENCURSOR;
  384.        INTRODUCTION;
  385.        NORMALCURSOR;
  386.      END;
  387.  
  388.    MAIN:
  389.     ClrScr;
  390.     FastCenter( 'Program SYNCSAT' , 2, LabelAttr );
  391.  
  392.    LOCATION:
  393.     FastWrite( 'Observer latitude ( -90° to +90° )     ' , 4, 5, TitleAttr );
  394.     gotoXY( 45, 4 );
  395.     Readln( OBSLAT );
  396.     OBSLATstr := Form( '###.##°', OBSLAT );
  397.  
  398.     FastWrite( 'Observer WEST Longitude ( 0° to 360° ) ' , 6,5, TitleAttr );
  399.     gotoXY( 45, 6 );
  400.     Readln( OBSLONG );
  401.     OBSLONGstr := Form( '###.##°', OBSLONG );
  402.  
  403.     FastWrite( 'Observer altitude ( meters ) ', 8,5, TitleAttr );
  404.     FastWrite( '< Above or below sea level > ', 9,5, HelpAttr );
  405.     gotoXY( 45, 8 );
  406.     Readln( OBSALT );
  407.     OBSALTstr := Form( '###.## m.', OBSALT );
  408.  
  409.                      { Read data from SAT.DAT file }
  410.    SATELLITE:
  411.     ASSIGN( FP, SATFILE );
  412.     RESET( FP );
  413.     J := 0;
  414.     WHILE NOT EOF( FP ) DO
  415.       Begin
  416.         WHILE NOT EOLN( FP ) DO
  417.           Begin
  418.             INC(J);
  419.             ReadLN( FP, SatRec[J].Bird, SatRec[J].Code, SatRec[J].GeoLong );
  420.           end; { while not eoln}
  421.       end; {while not EOF}
  422.     CLOSE( FP );
  423.     MaxRecs := J;
  424.  
  425.                   {Calculate results into record array }
  426.     For K := 1 to MaxRecs DO  { was 20 }
  427.       Begin
  428.         IF Str2Real( SatRec[K].Geolong, SATLONG ) THEN
  429.           Begin
  430.             Calculate_Obs_Position;
  431.             SatRec[K].Az :=  Form( '###.##°' , AZIMUTH );
  432.             SatRec[K].Az :=  SatRec[K].Az + AzQuad16( AZIMUTH );
  433.             SatRec[K].Alt:=  Form( '####.##°' , ELEVATION);
  434.           end;
  435.       end;
  436.  
  437.     HIDDENCURSOR;
  438.     Display_Results;
  439.     NORMALCURSOR;
  440.  
  441.                       { Request Another Selection }
  442.  
  443.     FastCenter( '                                 ', 25,    HeadAttr );
  444.     FastWrite( 'Another location  ( <Y>es, <N>o ) ' , 25,5, HelpAttr );
  445.     gotoXY( 40,25);
  446.     REPLY := ReadKey;
  447.     IF REPLY in [ 'Y', 'y' ] THEN
  448.      Begin
  449.        Clrscr;
  450.        goto LOCATION;
  451.      end
  452.     ELSE
  453.      EXIT;
  454.   end;
  455.  
  456.    { ********************* MAIN PROGRAM **********************************}
  457.   Begin
  458.     MainProgram;
  459.     ClrScr;
  460.   end.
  461.